### Out of Sight, Out of Mind? ###
### Figures in the Main Text ###

rm(list=ls())

library(magrittr) ## Version: 2.0.1
library(dplyr)  ## Version: 1.0.5
library(readstata13) ## Version: 0.9.2
library(interflex) ## Version: 1.2.5
library(ggplot2) ## Version: 3.3.3
library(gridExtra) ## Version: 2.3

d <- read.dta13('./final_data.dta')

## Figure 3: Proximity to and the congestion of primary health care services in Istanbul, before and after the Family Medicine Reform ----

plotdata <- d %>% filter(id_election==2 | id_election==4) %>% select(year, id_election, time_walk_min, pop1000_asm_dr)

## remove outliers for better visualization
plotdata_sub <- plotdata %>% filter(time_walk_min<30)

walkmin_density <- ggplot(plotdata_sub, aes(x=time_walk_min, fill = as.factor(year)))+geom_density(alpha = 0.5)+scale_fill_manual(name = 'Year', values=c("darkturquoise", "brown4"))+xlab('walking distance (in minutes)') + ylab('Density')+theme_bw()+theme(legend.position = 'bottom')

ggsave(walkmin_density, file = './figure3_left.pdf', width = 6, height = 5)


congestion_density <- ggplot(plotdata, aes(x=pop1000_asm_dr, fill = as.factor(year)))+geom_density(alpha = 0.5)+scale_fill_manual(name = 'Year', values=c("darkturquoise", "brown4"))+xlab('Number of people per doctor (in 1000s)') + ylab('Density')+theme_bw()+theme(legend.position = 'bottom')

ggsave(congestion_density, file = './figure3_right.pdf', width = 6, height = 5)


## Figure 4: Conditional marginal effects of DCongestion on DAKP using binning estimator ----

##(3 bins, the estimates are the medians of each bin). The moderator is property values in 2009 (logged). 

lograyic_models <- d %>% filter(year<2014) %>% select(school_id_cluster, Dakp_3, Dwalk_3, lograyic09, Duniversity_3, Dodr_3, Dydr_3, Dpopulation_3, Dhospital_pri_3, Dhospital_pub_3)

lograyic_models <- lograyic_models[complete.cases(lograyic_models), ]

walk_lograyic_binning <- interflex(Y = "Dakp_3", D = "Dwalk_3", X = "lograyic09", Z = c('Duniversity_3', 'Dodr_3', 'Dydr_3', 'Dpopulation_3', 'Dhospital_pri_3', 'Dhospital_pub_3'),data = lograyic_models,estimator = 'binning', cl="school_id_cluster", main = "", nbins = 3, Xlabel = "Property values (2009)", vcov.type = 'cluster', xlab = 'Moderator: Property values (2009)', ylab = 'Marginal effect of Dwalk on Dakp', cex.lab = 1, cex.axis = 1, theme.bw = TRUE, height=6, width = 4, file='./figure4.pdf', ylim = c(-0.35, 0.2))


## Figure 5: Conditional marginal effects of DCongestion on DAKP using binning estimator ----
##(3 bins, the estimates are the medians of each bin). The moderator is property values in 2009 (logged).

congestion_models <- d %>% filter(year<2014) %>% select(school_id_cluster, Dakp_3, Dpop1000_asm_dr_3, Dwalk_3, lograyic09, Duniversity_3, Dodr_3, Dydr_3, Dpopulation_3, Dhospital_pri_3, Dhospital_pub_3)

congestion_models <- congestion_models[complete.cases(congestion_models), ]


congestion_lograyic_binning <- interflex(Y = "Dakp_3", D = "Dpop1000_asm_dr_3", X = "lograyic09", Z = c('Duniversity_3', 'Dodr_3', 'Dydr_3', 'Dpopulation_3', 'Dhospital_pri_3', 'Dhospital_pub_3'),data = congestion_models,estimator = 'binning', cl="school_id_cluster", main = "", nbins = 3, Xlabel = "Property values (2009)", vcov.type = 'cluster', xlab = 'Moderator: Property values (2009)', ylab = 'Marginal effect of Dcongestion on Dakp', cex.lab = 1, cex.axis = 1, theme.bw = TRUE, height=6, width = 4, file='./figure5.pdf', ylim = c(-1.3, 0.4))


